Gravitational Collapse of Gravitational Waves in 3D Numerical Relativity 
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We demonstrate that evolutions of three-dimensional, strongly non-linear gravitational waves can 
be followed in numerical relativity, hence allowing many interesting studies of both fundamental 
and observational consequences. We study the evolution of time-symmetric, axisymmetric and non- 
axisymmetric Brill waves, including waves so strong that they collapse to form black holes under their 
own self-gravity. The critical amplitude for black hole formation is determined. The gravitational 
waves emitted in the black hole formation process are compared to those emitted in the head-on 
collision of two Misner black holes. 
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Gravitational waves have been an important area of 
research in Einstein's theory of gravity for years. Ein- 
stein's equations are nonlinear, and therefore can cause 
waves, which normally would disperse if weak enough, 
to be held together by their own gravity. This property 
characterizes Wheeler's geon |l],U proposed more than 
40 years ago, and is responsible for many interesting phe- 
nomena. Even in planar symmetric spacetimes, there are 
many interesting results, such as the formation of singu- 
larities from colliding plane waves (see [|| and references 
therein). In axisymmetry, Ref. studied the formation 
of black holes (BHs) by imploding gravitational waves, 
finding critical behavior B. 

These discoveries are all in spacetimes with special 
symmetries, but they raise important questions about 
fully general three-dimensional (3D) spacetimes, e.g. the 
nature of critical phenomena in the absence of symme- 
tries has only recently been studied through a perturba- 
tive approach 3D studies of fully nonlinear gravity 
can only be made through the machinery of numerical 
relativity. A few studies of gravitational wave evolu- 
tions have been performed in the linear and near linear 
regimes j7j-|J, in preparation for the study of fully non- 
linear, strong field 3D wave dynamics. However, until 
now no such studies have been successfully carried out. 

In this paper we present the first successful simula- 
tions of highly nonlinear gravitational waves in 3D. We 
study the process of strong waves collapsing to form BHs 
under their own self-gravity We determine the critical 
amplitude for the formation of BHs and show that one 
can now carry out these evolutions for long times. For 
waves that are not strong enough to form BHs, we follow 
their implosion, bounce and dispersal. For waves strong 
enough to collapse to a BH under their own self gravity, 
we find the dynamically formed apparent horizons (AHs), 
and extract the gravitational radiation generated in the 
collapse process. These waveforms can be compared in 



axisymmetry to head-on BH collisions (performed earlier 
and reported in ]To[|). The waveforms are similar at late 
times, dominated by the quasi-normal modes of the re- 
sulting BHs as expected. The difference in the waveforms 
at early times for these two very different collapse sce- 
narios shows to what extent one can extract information 
about the BH formation process from the observation of 
the gravitational radiation emitted by the system. All 
the simulations presented here were performed with the 
newly developed Cactus code. For a description of the 
code and the numerical methods used, s ee [|ll| - |l5| . 

We take as initial data a pure Brill fll6fl ty pe gravi- 
tational wave, later studied by Eppley 117 18 and oth- 
ers ilStl. The metric takes the form 



ds 2 = * 4 [e 2 ? {dp 2 + dz 2 ) + p 2 d<f> 2 ] = *W 



(1) 



where q is a free function subject to certain boundary 
conditions. Following p4] , p0 21 1, we choose q of the form 
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where a, c are constants, r 2 = p 2 + z 2 and n is an inte- 
ger. For c = 0, these data sets reduce to the Holz M 
axisymmetric form, recently studied in full 3D Carte- 
sian coordinates in preparation for the present work [p2[ . 
Taking this form for q, we impose the condition of time- 
symmetry, and solve the Hamiltonian constraint numeri- 
cally in Cartesian coordinates. An initial data set is thus 
characterized only by the parameters (a, c, n). For the 
case (a, 0, 0), we found in that no AH exists in initial 
data for a < 11.8, and we also studied the appearance of 
an AH for other values of c and n. 

Such initial data can be evolved in full 3D using the 
Cactus code, which allows the use of different formula- 
tions of the Einstein equations, different coordinate con- 
ditions, and different numerical methods. Our focus here 
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is on new physics, but since stable evolutions of such 
strong gravitational waves have not been obtained be- 
fore, we comment briefly on the method used for the 
results in this paper. In ply, Baumgarte and Shapiro 
note for weak waves that a system, which is essentially 
the conformally decomposed ADM system of Shibata and 
Nakamura shows greatly increased numerical stabil- 
ity over the standard ADM formulation p4| ]. We will 
refer to this system as the BSSN formulation. The use of 
a particular connection variable in BSSN is reminiscent 
of the Bona-Masso formulation [^5 26 1. We found that 
BSSN as given in [p3| with maximal slicing, a 3-step it- 
erative Crank-Nicholson (ICN) scheme, and a radiative 
(Sommerfeld) boundary condition is very stable and reli- 
able even for the strong waves considered here. The key 
new extensions to previous BSSN results are that the 
stability can be extended to (i) strong, dynamical fields 
and (ii) maximal slicing, where the latter requires some 
care. Maximal slicing is defined by vanishing of the mean 
extrinsic curvature, K=0, and the BSSN formulation al- 
lowed us to cleanly implement this feature numerically, 
in contrast with the standard ADM equations. (A re- 
lated idea to improve stability with maximal slicing is 
that of K-drivers, which helps dramatically, but is ulti- 
mately not sufficient for very strong waves in standard 
ADM formulations [^tJ, but compare [p6|.) 

We begin our discussion of the physical results with 
the parameter set (a=4, c—0, n=0); a rather strong ax- 
isymmetric Brill wave (BW). Even though this data set is 
axisymmctric, the evolution has been carried out in full 
3D, exploiting the reflection symmetry on the coordinate 
planes to evolve only one of the eight octants. The evo- 
lution of this data set shows that part of the wave prop- 
agates outward while part implodes, re-expanding after 
passing through the origin. However, due to the non- 
linear self-gravity, not all of it immediately disperses out 
to infinity; again part re-collapses and bounces again. Af- 
ter a few collapses and bounces the wave completely dis- 
perses out to infinity. This behavior is shown in Fig. ^a, 
where the evolution of the central value of the lapse 
is given for simulations with three different grid sizes: 
Ax—Ay=Az=0.16 (low resolution), 0.08 (medium res- 
olution) and 0.04 (high resolution), using 32 3 , 64 3 and 
128 3 grid points respectively. At late times, the lapse 
returns to 1 (the log returns to 0). Fig. |l]b shows the 
evolution of the log of the central value of the Riemann 
invariant J for the same resolutions. At late times J set- 
tles on a constant value that converges rapidly to zero as 
we refine the grid. With these results, and direct veri- 
fication that the metric functions become stationary at 
late times, we conclude that spacetime returns to flat (in 
non-trivial spatial coordinates; the metric is decidedly 
non-flat in appearance!). 

Next we increase the amplitude to a = 6, holding other 
parameters fixed. Fig. || shows the evolution of the lapse 
and the Riemann invariant J for this case, showing a 
clear contrast with Fig. |l|. The lapse now collapses imme- 
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FIG. 1. a) Evolution of the log of the lapse a at r=0 
for the axisymmetric data (4,0,0). The dashed/dotted/solid 
lines represent simulations at low/medium/high resolution, 
b) Evolution of the Riemann invariant J at r=0. The wave 
disperses after dynamic evolution, leaving flat space behind. 



diately, and the Riemann invariant after an initial drop 
grows to a large value at the origin until it is halted by the 
collapse of the lapse. For this amplitude the low resolu- 
tion is now too crude and the code crashes at t ~ 10. We 
have therefore added an extra simulation with Aa;=0.053 
("intermediate" resolution) using 96 3 grid points. 

To confirm that a BH has indeed formed, we searched 
for an AH in the a=6 case (using a minimization algo- 
rithm P2| ). For high resolution, an AH was first found 
at t=7.7, which grows slowly in both coordinate radius 
and area. Fig. || shows the location of the AH on the 
x-z plane at time t = 10 for the three resolutions. The 
mass of the horizon at this time is about Mah = 0.87, 
but then due to poor resolution of the grid stretching 
(a common problem of all BH simulations with singu- 
larity avoiding slicings), it continues to grow, ultimately 
exceeding the initial ADM mass of the spacetime, which 
for this data set is Madm = 0.99 (obtained in the way 
described in p2||). However, the total energy radiated 
is about 0.12, computed from the Zerilli functions, com- 
pletely consistent with Mah = 0.87 and an initial mass 
of Madm = 0.99 . CPU time constraints make it difficult 
to run long term, higher resolution simulations (high res- 
olution used ~ 120 hours running on 16 processors of an 
SGI/Cray-Origin 2000). We also confirmed that an event 
horizon does not exist in the initial data by integrating 
null surfaces out from the origin during the simulation. 

From these two studies we conclude that the critical 
amplitude a* for BH formation for the axisymmetric BW 
packet is a* = 5± 1. We have performed more simu- 
lations within this range, and have narrowed down the 
interval to a* = 4.85 ±0.15, although near the critical 
solution higher resolution is required to establish conver- 
gence. Our study of these near-critical solutions is still 
under way and will be presented elsewhere. 

It is particularly exciting that the dynamical evolution 
can be followed long enough for the extraction of gravi- 
tational waveforms even for the BH formation case. One 
important question is what physical information of the 
gravitational collapse process can be extracted from the 
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FIG. 2. a) Evolution of the lapse a at r=0 for the ax- 
isymmetric data set (6,0,0). The dashed, dotted, solid and 
dashed-dotted lines represent simulations at low, medium, in- 
termediate and high resolutions respectively, b) Evolution of 
the Riemann invariant J at r=0. c) Coordinate location of 
the dynamically formed AH on the x-z plane at t=10. 



observation of the radiation. How much will the wave- 
forms from different BH formation processes be different? 
For this purpose we compare the BW collapse waveforms 
to those of a very different collapse process, namely the 
head-on collision of two BHs. In Fig. || we show the 
{l=2,m—0} Zcrilli function ip, obtained from the evolu- 
tion of Misner data for n = 1.2,1.8,2.2 JTo) , and from 
the axisymmetric a — 6 BW collapse. (The case ^ = 1.2 
represents a single perturbed black hole, at ^ = 2.2 there 
are two separate black holes that are outside the pertur- 
bative regime.) To compare the waveforms, we adjust the 
time coordinate of the BW waveforms based on the time 
delay for different "detector" positions, which for the BW 
is at r = \ AM adm and for the BHs at r — 20Madm- 
We also scale the Zerilli function amplitude for the BHs 
by Madm and the BW by IOMadm to put them on the 
same figure. 

We notice the following: (1) The BW waveform is dom- 
inated by quasi-normal modes (QNMs) at late times just 
like in the 2BH case, as expected. A QNM fit shows 
that at about lQM a d m from the beginning of the wave- 
train the fundamental mode dominates. (2) However, 
the BW waveform has more high frequency QNM com- 
ponents in the early phase. The waveforms start with a 
different offset from zero, which is substantially larger in 
magnitude in the Brill wave case, but note that in the 
BW case the detector is put much closer in (at 4.6M a£ ; m ) 
and the Zerilli function extraction process gives a 

larger "Coulomb" component [29[| . (3) The fundamental 
QNMs that dominate the late time evolutions for the two 
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FIG. 3. We compare the l=2,m—0 extracted waveform 
for the head-on collision of BHs obtained by |l(]] of the 
fj, — 1.2,1.8,2.2 Misner data (solid lines, increasing ampli- 
tude) to that of the a — 6 collapsing BW (dotted and dashed 
lines). 



cases have the same phase! We see that both waveforms 
dip at 30M a d m , and peak at 38M a d m , to high accuracy. 
We note that the 2BH waveforms for all /i < 2.2 have 
their fundamental QNM appearing with about the same 
phase, and we see here the BW collapse case also has the 
same phase. This and other interesting comparisons be- 
tween the two collapse scenarios will be discussed further 
elsewhere. (We note that these features noted above are 
not sensitive to the fi value chosen, within the range of 
fx = 1.2 - 2.2.) 

Next we go to a pure strong wave case with full 3D 
features (the first ever simulated) , where the initial wave- 
form is even more dominated by details of the BH forma- 
tion process. Fig. [| shows the development of the data 
set (<z=6, c=0.2, n—1), which has reflection symmetry 
across coordinate planes; it again suffices to evolve only 
an octant. The initial ADM mass of this data set turns 
out to be Madm = 1.12. Fig. ||a shows a comparison of 
the AHs of this 3D and the previous axisymmetric cases, 
using the same high resolution, at t=10 on the x-z plane. 
The mass of the 3D AH case is larger, weighing in at 
M AH =Q.m (compared to M AH (2D) = 0.87). 

In Fig. |]b we show the {/=2,m=0} waveform of this 3D 
case, compared to the previous axisymmetric case. The 
c = 0.2 waveform has a longer wave length at late times, 
consistent with the fact that a larger mass BH is formed 
in the 3D case. Figs. ||c and||d show the same comparison 
for the {/=4,m=0} and {7=2, m=2} modes respectively. 
Notice that while the first two modes are of similar am- 
plitude for both runs, the 3D {7=2, m=2} mode is com- 
pletely different; as a non-axisymmetric contribution, it is 
absent in the axisymmetric run (in fact, it doesn't quite 
vanish due to numerical error, but it remains of order 
10~ 6 ). We also show a fit to the corresponding QNM's 
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of a BH of mass 1.0. The fit was performed in the time 
interval (10,36), and is noticeably worse if the fit is at- 
tempted to earlier times, again showing that the lowest 
QNM's dominate at around 10. The early parts of the 
waveforms t < 10 reflect the details of the initial data 
and BH formation process. This is especially clear in 
the {/=2,m=2} mode, which seems to provide the most 
information about the initial data and the 3D BH forma- 
tion process. At present no 3D BH formation simulation 
from other scenarios (e.g., true spiraling BH coalescence) 
are available for comparison, as in the axisymmetric case, 
but such simulations may actually be available soon |j0[ . 
It will be interesting to compare such studies with 3D 
wave collapses, such as that presented here. 

In conclusion, we demonstrated numerical evolutions 
of 3D, strongly non-linear gravitational waves, and stud- 
ied gravitational collapse of axisymmetric and non- 
axisymmctric gravitational waves. We compared the 
wave collapse to the head-on collision of two black holes. 
The research opens the door to many investigations. 
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FIG. 4. a) The solid (dotted) line is the AH for the full 
3D data set (6, 0.2, 1) ((6, 0, 0)) at t=10 on the x-z plane, b) 
The {Z=2,m=0} waveform for the 3D (6, 0.2, 1) case at r = 4 
(solid line) is compared to axisymmetric (6, 0, 0) case (dotted 
line) . The dashed line shows the fit of the 3D case to the cor- 
responding mode for a BH of mass 1.0. c) Same comparison 
for the {7=4, m=0} waveform, d) Same comparison for the 
non-axisymmetric {Z=2,m=2} waveform. 
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